df_ger_covid <- read_csv('Germany_timeseries_prep.csv')
df_ger_socdist <- read_csv('Germany_socdist_fb_kreis.csv')
df_ger_ctrl <- read_delim('Germany_controls.csv', delim = ';')
# prevalence
df_ger_covid_clean <- df_ger_covid %>% mutate(date = as.Date(date, "%d%b%Y"),
kreis = as.character(kreis)) %>%
filter(date >= '2020-03-09' & date <= '2020-03-31') %>%
group_by(kreis) %>%
mutate(time = row_number()) %>%
ungroup() %>%
dplyr::select(-runday, -kreis_name, -ewz, -shape__area,
-cumcase, -anzahlfall, -popdens) %>%
dplyr::rename(pers_o = open,
pers_c = sci,
pers_e = extra,
pers_a = agree,
pers_n = neuro)
df_ger_covid_clean %>% head()
# social distancing
df_ger_socdist_clean <- df_ger_socdist %>%
filter(date >= '2020-03-09' & date <= '2020-03-31') %>%
mutate(kreis = as.character(kreis)) %>%
group_by(kreis) %>%
mutate(time = row_number()) %>%
ungroup() %>%
rename(socdist_single_tile = all_day_ratio_single_tile_users) %>%
select(kreis, time, socdist_single_tile)
df_ger_socdist_clean %>% head()
NA
# controls
df_ger_ctrl_clean <- df_ger_ctrl %>% select(-kreis_nme) %>%
mutate(kreis = as.character(kreis),
hospital_beds = as.numeric(str_replace(hospital_beds, ',', '.')))
df_ger_ctrl_clean %>% head()
NA
NA
# merge
df_ger <- df_ger_covid_clean %>%
plyr::join(df_ger_socdist_clean, by = c('kreis', 'time')) %>%
inner_join(df_ger_ctrl_clean, by = 'kreis')
df_ger
NA
df_ger %>% ggplot(aes(x=time, y=rate_day)) +
geom_point(aes(col=kreis, size=popdens)) +
geom_smooth(method="loess", se=T) +
theme(legend.position="none") +
ggtitle("Overall prevalence over time")
pers <- c('pers_o', 'pers_c', 'pers_e', 'pers_a', 'pers_n')
for (i in pers){
gg <- df_ger %>% mutate(prev_tail = cut(.[[i]],
breaks = c(-Inf, quantile(.[[i]], 0.2), quantile(.[[i]], 0.8), Inf),
labels = c('lower tail', 'center', 'upper tail'))) %>%
filter(prev_tail != 'center') %>%
ggplot(aes(x=time, y=rate_day)) +
geom_point(aes(col=kreis, size=popdens)) +
geom_smooth(method="loess", se=T) +
facet_wrap(~prev_tail) +
theme(legend.position="none") +
ggtitle(i)
print(gg)
}
weekend <- c(6, 7, 13, 14, 20, 21)
df_ger_loess <- df_ger %>% filter(!time %in% weekend) %>%
split(.$kreis) %>%
map(~ loess(socdist_single_tile ~ time, data = .)) %>%
map(predict, 1:23) %>%
bind_rows() %>%
gather(key = 'kreis', value = 'loess') %>%
group_by(kreis) %>%
mutate(time = row_number())
df_ger <- df_ger %>% merge(df_ger_loess, by=c('kreis', 'time')) %>%
mutate(socdist_single_tile_clean = ifelse(time %in% weekend, loess,
socdist_single_tile)) %>%
arrange(kreis, time)
df_ger %>% ggplot(aes(x=time, y=loess, group=kreis)) +
geom_line()
NA
NA
df_ger %>% ggplot(aes(x=time, y=socdist_single_tile_clean)) +
geom_point(aes(col=kreis, size=popdens)) +
geom_smooth(method="loess", se=T) +
theme(legend.position="none") +
ggtitle("Overall social distancing (single tile) over time")
pers <- c('pers_o', 'pers_c', 'pers_e', 'pers_a', 'pers_n')
for (i in pers){
gg <- df_ger %>% mutate(prev_tail = cut(.[[i]],
breaks = c(-Inf, quantile(.[[i]], 0.2), quantile(.[[i]], 0.8), Inf),
labels = c('lower tail', 'center', 'upper tail'))) %>%
filter(prev_tail != 'center') %>%
ggplot(aes(x=time, y=socdist_single_tile_clean)) +
geom_point(aes(col=kreis, size=popdens)) +
geom_smooth(method="loess", se=T) +
facet_wrap(~prev_tail) +
theme(legend.position="none") +
ggtitle(i)
print(gg)
}
df_ger %>% group_by(kreis) %>%
summarize_if(is.numeric, mean) %>%
select(-kreis, -time) %>%
cor(use='pairwise.complete') %>%
round(3)
pers_e pers_a pers_c pers_n pers_o rate_day socdist_single_tile women academics cdu afd hospital_beds
pers_e 1.000 0.223 0.255 -0.375 0.277 0.168 0.149 -0.011 0.162 0.111 -0.108 -0.037
pers_a 0.223 1.000 0.347 -0.379 0.167 0.135 0.002 -0.009 0.180 0.023 0.015 0.004
pers_c 0.255 0.347 1.000 -0.373 -0.063 -0.010 0.030 -0.046 -0.045 -0.006 0.140 -0.157
pers_n -0.375 -0.379 -0.373 1.000 -0.046 -0.127 -0.141 -0.012 -0.068 -0.061 0.064 0.111
pers_o 0.277 0.167 -0.063 -0.046 1.000 0.115 0.090 -0.094 0.431 -0.053 -0.187 0.281
rate_day 0.168 0.135 -0.010 -0.127 0.115 1.000 0.265 -0.026 0.128 0.136 -0.135 -0.060
socdist_single_tile 0.149 0.002 0.030 -0.141 0.090 0.265 1.000 0.049 -0.018 0.149 -0.301 -0.243
women -0.011 -0.009 -0.046 -0.012 -0.094 -0.026 0.049 1.000 -0.001 0.074 -0.089 0.009
academics 0.162 0.180 -0.045 -0.068 0.431 0.128 -0.018 -0.001 1.000 -0.208 -0.117 0.335
cdu 0.111 0.023 -0.006 -0.061 -0.053 0.136 0.149 0.074 -0.208 1.000 -0.157 -0.091
afd -0.108 0.015 0.140 0.064 -0.187 -0.135 -0.301 -0.089 -0.117 -0.157 1.000 -0.030
hospital_beds -0.037 0.004 -0.157 0.111 0.281 -0.060 -0.243 0.009 0.335 -0.091 -0.030 1.000
tourism_beds -0.115 -0.134 -0.102 0.007 -0.108 -0.038 -0.012 0.003 -0.156 0.115 -0.011 -0.043
gdp 0.121 0.064 0.000 -0.067 0.247 0.110 -0.047 -0.044 0.374 -0.078 -0.133 0.266
manufact -0.009 -0.061 -0.051 -0.007 -0.038 0.076 -0.071 0.026 -0.088 0.019 0.065 0.070
airport -0.188 -0.167 -0.112 0.171 -0.222 -0.031 -0.142 0.043 -0.331 0.105 0.183 -0.046
age 0.007 0.086 0.008 -0.016 -0.063 -0.034 -0.017 -0.046 -0.051 0.020 0.091 0.029
popdens 0.062 0.038 0.018 -0.080 0.072 0.092 0.110 -0.001 0.056 -0.016 -0.028 0.024
loess 0.139 -0.013 0.020 -0.143 0.105 0.229 0.979 0.040 0.012 0.124 -0.344 -0.229
socdist_single_tile_clean 0.138 -0.015 0.020 -0.142 0.104 0.228 0.979 0.041 0.010 0.123 -0.344 -0.228
tourism_beds gdp manufact airport age popdens loess socdist_single_tile_clean
pers_e -0.115 0.121 -0.009 -0.188 0.007 0.062 0.139 0.138
pers_a -0.134 0.064 -0.061 -0.167 0.086 0.038 -0.013 -0.015
pers_c -0.102 0.000 -0.051 -0.112 0.008 0.018 0.020 0.020
pers_n 0.007 -0.067 -0.007 0.171 -0.016 -0.080 -0.143 -0.142
pers_o -0.108 0.247 -0.038 -0.222 -0.063 0.072 0.105 0.104
rate_day -0.038 0.110 0.076 -0.031 -0.034 0.092 0.229 0.228
socdist_single_tile -0.012 -0.047 -0.071 -0.142 -0.017 0.110 0.979 0.979
women 0.003 -0.044 0.026 0.043 -0.046 -0.001 0.040 0.041
academics -0.156 0.374 -0.088 -0.331 -0.051 0.056 0.012 0.010
cdu 0.115 -0.078 0.019 0.105 0.020 -0.016 0.124 0.123
afd -0.011 -0.133 0.065 0.183 0.091 -0.028 -0.344 -0.344
hospital_beds -0.043 0.266 0.070 -0.046 0.029 0.024 -0.229 -0.228
tourism_beds 1.000 -0.094 -0.051 0.297 0.075 -0.157 0.031 0.031
gdp -0.094 1.000 0.339 -0.085 -0.014 0.082 -0.026 -0.027
manufact -0.051 0.339 1.000 0.136 -0.049 0.072 -0.119 -0.120
airport 0.297 -0.085 0.136 1.000 0.077 -0.179 -0.169 -0.169
age 0.075 -0.014 -0.049 0.077 1.000 -0.038 -0.021 -0.021
popdens -0.157 0.082 0.072 -0.179 -0.038 1.000 0.094 0.094
loess 0.031 -0.026 -0.119 -0.169 -0.021 0.094 1.000 1.000
socdist_single_tile_clean 0.031 -0.027 -0.120 -0.169 -0.021 0.094 1.000 1.000
# function calculates all relevant models
run_models <- function(y, lvl1_x, lvl2_x, lvl2_id, data, ctrls=F){
# subset data
data = data %>%
dplyr::select(all_of(y), all_of(lvl1_x), all_of(lvl2_x), all_of(lvl2_id),
popdens, rate_day)
data = data %>%
dplyr::rename(y = all_of(y),
lvl1_x = all_of(lvl1_x),
lvl2_x = all_of(lvl2_x),
lvl2_id = all_of(lvl2_id)
)
# configure optimization procedure
ctrl_config <- lmeControl(opt = 'optim', maxIter = 100, msMaxIter = 100)
# baseline
baseline <- lme(fixed = y ~ 1, random = ~ 1 | lvl2_id,
data = data,
correlation = corAR1(),
control = ctrl_config,
method = 'ML')
# random intercept fixed slope
random_intercept <- lme(fixed = y ~ lvl1_x + lvl2_x,
random = ~ 1 | lvl2_id,
data = data,
correlation = corAR1(),
control = ctrl_config,
method = 'ML')
# random intercept random slope
random_slope <- lme(fixed = y ~ lvl1_x + lvl2_x,
random = ~ lvl1_x | lvl2_id,
data = data,
correlation = corAR1(),
control = ctrl_config,
method = 'ML')
# cross level interaction
interaction <- lme(fixed = y ~ lvl1_x * lvl2_x,
random = ~ lvl1_x | lvl2_id,
data = data,
correlation = corAR1(),
control = ctrl_config,
method = 'ML')
# create list with results
results <- list('baseline' = baseline,
"random_intercept" = random_intercept,
"random_slope" = random_slope,
"interaction" = interaction)
if (ctrls == 'dem' | ctrls == 'prev'){
# random intercept random slope
random_slope_ctrl_dem <- lme(fixed = y ~ lvl1_x + lvl2_x + popdens,
random = ~ lvl1_x | lvl2_id,
data = data,
correlation = corAR1(),
control = ctrl_config,
method = 'ML')
# cross level interaction
interaction_ctrl_main_dem <- lme(fixed = y ~ lvl1_x * lvl2_x + popdens,
random = ~ lvl1_x | lvl2_id,
data = data,
correlation = corAR1(),
control = ctrl_config,
method = 'ML')
# cross level interaction
interaction_ctrl_int_dem <- lme(fixed = y ~ lvl1_x * lvl2_x + lvl1_x * popdens,
random = ~ lvl1_x | lvl2_id,
data = data,
correlation = corAR1(),
control = ctrl_config,
method = 'ML')
# create list with results
results <- list('baseline' = baseline,
"random_intercept" = random_intercept,
"random_slope" = random_slope,
"interaction" = interaction,
"random_slope_ctrl_dem" = random_slope_ctrl_dem,
"interaction_ctrl_main_dem" = interaction_ctrl_main_dem,
"interaction_ctrl_int_dem" = interaction_ctrl_int_dem)
}
if (ctrls == 'prev'){
# random intercept random slope
random_slope_ctrl_prev <- lme(fixed = y ~ lvl1_x + lvl2_x + popdens + rate_day,
random = ~ lvl1_x + rate_day | lvl2_id,
data = data,
correlation = corAR1(),
control = ctrl_config,
method = 'ML')
# cross level interaction
interaction_ctrl_main_prev <- lme(fixed = y ~ lvl1_x * lvl2_x + popdens + rate_day,
random = ~ lvl1_x | lvl2_id,
data = data,
correlation = corAR1(),
control = ctrl_config,
method = 'ML')
# cross level interaction
interaction_ctrl_int_prev<- lme(fixed = y ~ lvl1_x * lvl2_x + lvl1_x * popdens + rate_day,
random = ~ lvl1_x + rate_day | lvl2_id,
data = data,
correlation = corAR1(),
control = ctrl_config,
method = 'ML')
# create list with results
results <- list('baseline' = baseline,
"random_intercept" = random_intercept,
"random_slope" = random_slope,
"interaction" = interaction,
"random_slope_ctrl_dem" = random_slope_ctrl_dem,
"interaction_ctrl_main_dem" = interaction_ctrl_main_dem,
"interaction_ctrl_int_dem" = interaction_ctrl_int_dem,
"random_slope_ctrl_prev" = random_slope_ctrl_prev,
"interaction_ctrl_main_prev" = interaction_ctrl_main_prev,
"interaction_ctrl_int_prev" = interaction_ctrl_int_prev)
}
if(ctrls == 'exp'){
# random intercept random slope
random_slope_exp <- lme(fixed = y ~ (lvl1_x + I(lvl1_x^2)) + lvl2_x,
random = ~ (lvl1_x + I(lvl1_x^2)) | lvl2_id,
data = data,
correlation = corAR1(),
control = ctrl_config,
method = 'ML')
# cross level interaction
interaction_exp <- lme(fixed = y ~ (lvl1_x + I(lvl1_x^2)) * lvl2_x,
random = ~ (lvl1_x + I(lvl1_x^2)) | lvl2_id,
data = data,
correlation = corAR1(),
control = ctrl_config,
method = 'ML')
# create list with results
results <- list('baseline' = baseline,
"random_intercept" = random_intercept,
"random_slope" = random_slope,
"interaction" = interaction,
"random_slope_exp" = random_slope_exp,
"interaction_exp" = interaction_exp)
}
return(results)
}
# extracts table with coefficients and tests statistics
extract_results <- function(models) {
models_summary <- models %>%
map(summary) %>%
map("tTable") %>%
map(as.data.frame) %>%
map(round, 10)
# %>% map(~ .[str_detect(rownames(.), 'Inter|lvl'),])
return(models_summary)
}
compare_models <- function(models) {
mdl_names <- models %>% names()
str = ''
for (i in mdl_names){
mdl_str <- paste('models$', i, sep = '')
if(i == 'baseline'){
str <- mdl_str
}else{
str <- paste(str, mdl_str, sep=', ')
}
}
anova_str <- paste0('anova(', str, ')')
mdl_comp <- eval(parse(text=anova_str))
rownames(mdl_comp) = mdl_names
return(mdl_comp)
}
df_ger_scaled <- df_ger %>% dplyr::select(kreis, time, pers_o,
pers_c, pers_e, pers_e, pers_a, pers_n,
popdens, rate_day, socdist_single_tile) %>%
mutate_at(vars(-kreis, -time), scale)
df_ger_scaled %>% head()
models_o_covid <-run_models(y = 'rate_day',
lvl1_x = 'time',
lvl2_x = 'pers_o',
lvl2_id = 'kreis',
data = df_ger_scaled,
ctrls = 'dem')
extract_results(models_o_covid)
$baseline
$random_intercept
$random_slope
$interaction
$random_slope_ctrl_dem
$interaction_ctrl_main_dem
$interaction_ctrl_int_dem
compare_models(models_o_covid)
NA
models_c_covid <-run_models(y = 'rate_day',
lvl1_x = 'time',
lvl2_x = 'pers_c',
lvl2_id = 'kreis',
data = df_ger_scaled,
ctrls = 'dem')
extract_results(models_c_covid)
$baseline
$random_intercept
$random_slope
$interaction
$random_slope_ctrl_dem
$interaction_ctrl_main_dem
$interaction_ctrl_int_dem
compare_models(models_c_covid)
NA
NA
models_e_covid <-run_models(y = 'rate_day',
lvl1_x = 'time',
lvl2_x = 'pers_e',
lvl2_id = 'kreis',
data = df_ger_scaled,
ctrls = 'dem')
extract_results(models_e_covid)
$baseline
$random_intercept
$random_slope
$interaction
$random_slope_ctrl_dem
$interaction_ctrl_main_dem
$interaction_ctrl_int_dem
compare_models(models_e_covid)
NA
NA
models_a_covid <-run_models(y = 'rate_day',
lvl1_x = 'time',
lvl2_x = 'pers_a',
lvl2_id = 'kreis',
data = df_ger_scaled,
ctrls = 'dem')
extract_results(models_a_covid)
$baseline
$random_intercept
$random_slope
$interaction
$random_slope_ctrl_dem
$interaction_ctrl_main_dem
$interaction_ctrl_int_dem
compare_models(models_a_covid)
NA
NA
models_n_covid <-run_models(y = 'rate_day',
lvl1_x = 'time',
lvl2_x = 'pers_n',
lvl2_id = 'kreis',
data = df_ger_scaled,
ctrls = 'dem')
extract_results(models_n_covid)
$baseline
$random_intercept
$random_slope
$interaction
$random_slope_ctrl_dem
$interaction_ctrl_main_dem
$interaction_ctrl_int_dem
compare_models(models_n_covid)
NA
NA
models_o_covid_exp <-run_models(y = 'rate_day',
lvl1_x = 'time',
lvl2_x = 'pers_o',
lvl2_id = 'kreis',
data = df_ger_scaled,
ctrls = 'exp')
extract_results(models_o_covid_exp)
$baseline
$random_intercept
$random_slope
$interaction
$random_slope_exp
$interaction_exp
compare_models(models_o_covid_exp)
NA
models_c_covid_exp <-run_models(y = 'rate_day',
lvl1_x = 'time',
lvl2_x = 'pers_c',
lvl2_id = 'kreis',
data = df_ger_scaled,
ctrls = 'exp')
extract_results(models_c_covid_exp)
$baseline
$random_intercept
$random_slope
$interaction
$random_slope_exp
$interaction_exp
compare_models(models_c_covid_exp)
NA
models_e_covid_exp <-run_models(y = 'rate_day',
lvl1_x = 'time',
lvl2_x = 'pers_e',
lvl2_id = 'kreis',
data = df_ger_scaled,
ctrls = 'exp')
extract_results(models_e_covid_exp)
$baseline
$random_intercept
$random_slope
$interaction
$random_slope_exp
$interaction_exp
compare_models(models_e_covid_exp)
NA
models_a_covid_exp <-run_models(y = 'rate_day',
lvl1_x = 'time',
lvl2_x = 'pers_a',
lvl2_id = 'kreis',
data = df_ger_scaled,
ctrls = 'exp')
extract_results(models_a_covid_exp)
$baseline
$random_intercept
$random_slope
$interaction
$random_slope_exp
$interaction_exp
compare_models(models_a_covid_exp)
NA
models_n_covid_exp <-run_models(y = 'rate_day',
lvl1_x = 'time',
lvl2_x = 'pers_n',
lvl2_id = 'kreis',
data = df_ger_scaled,
ctrls = 'exp')
extract_results(models_n_covid_exp)
$baseline
$random_intercept
$random_slope
$interaction
$random_slope_exp
$interaction_exp
compare_models(models_n_covid_exp)
NA
summary_table <- function(models, dv_name){
temp_df_ctrl_main <- NULL
temp_df_ctrl_int <- NULL
for (i in models){
results <- i %>% extract_results()
results_ctrl_main <- results$interaction_ctrl_main_dem['lvl1_x:lvl2_x',]
temp_df_ctrl_main <- temp_df_ctrl_main %>% rbind(results_ctrl_main)
results_ctrl_int <- results$interaction_ctrl_int_dem['lvl1_x:lvl2_x',]
temp_df_ctrl_int <- temp_df_ctrl_int %>% rbind(results_ctrl_int)
}
names_ctrl_main <- paste0(dv_name, '~', c('o', 'c', 'e', 'a', 'n'), '*time', '_crtl_popdens')
names_ctrl_int <- paste0(dv_name, '~', c('o', 'c', 'e', 'a', 'n'), '*time', '_crtl_popdens*time')
rownames(temp_df_ctrl_main) <- names_ctrl_main
rownames(temp_df_ctrl_int) <- names_ctrl_int
sum_tab <- rbind(temp_df_ctrl_main, temp_df_ctrl_int) %>% round(4)
return(sum_tab)
}
# prevalence
models_prev <- list(models_o_covid,
models_c_covid,
models_e_covid,
models_a_covid,
models_n_covid)
sum_tab_prev <- summary_table(models_prev, dv_name = 'prev')
write.table(sum_tab_prev, quote=F)
Value Std.Error DF t-value p-value
prev~o*time_crtl_popdens 0.0071 0.0029 8798 2.4364 0.0149
prev~c*time_crtl_popdens 3e-04 0.003 8798 0.0944 0.9248
prev~e*time_crtl_popdens 0.0099 0.0029 8798 3.402 7e-04
prev~a*time_crtl_popdens 0.0091 0.0029 8798 3.1314 0.0017
prev~n*time_crtl_popdens -0.0072 0.0029 8798 -2.4642 0.0138
prev~o*time_crtl_popdens*time 0.0069 0.0029 8797 2.3439 0.0191
prev~c*time_crtl_popdens*time 2e-04 0.0029 8797 0.0687 0.9452
prev~e*time_crtl_popdens*time 0.0097 0.0029 8797 3.3233 9e-04
prev~a*time_crtl_popdens*time 0.009 0.0029 8797 3.086 0.002
prev~n*time_crtl_popdens*time -0.0069 0.0029 8797 -2.3608 0.0183
# social distancing
models_socdist <- list(models_o_socdist,
models_c_socdist,
models_e_socdist,
models_a_socdist,
models_n_socdist)
sum_tab_socdist <- summary_table(models_socdist, dv_name = 'socdist')
write.table(sum_tab_socdist, quote=F)
Value Std.Error DF t-value p-value
socdist~o*time_crtl_popdens 0.0041 0.0014 8798 3.0071 0.0026
socdist~c*time_crtl_popdens 0 0.0014 8798 -0.0081 0.9935
socdist~e*time_crtl_popdens 0.003 0.0014 8798 2.1539 0.0313
socdist~a*time_crtl_popdens 0.0017 0.0014 8798 1.2497 0.2114
socdist~n*time_crtl_popdens -0.0015 0.0014 8798 -1.0885 0.2764
socdist~o*time_crtl_popdens*time 0.0041 0.0014 8797 2.9445 0.0032
socdist~c*time_crtl_popdens*time 0 0.0014 8797 -0.0257 0.9795
socdist~e*time_crtl_popdens*time 0.0029 0.0014 8797 2.0978 0.036
socdist~a*time_crtl_popdens*time 0.0017 0.0014 8797 1.2133 0.2251
socdist~n*time_crtl_popdens*time -0.0014 0.0014 8797 -1.0139 0.3106
# slope prevalence
df_ger_slope_prev <- df_ger %>% split(.$kreis) %>%
map(~ lm(rate_day ~ time, data = .)) %>%
map(coef) %>%
map_dbl('time') %>%
as.data.frame() %>%
rownames_to_column('kreis') %>%
rename(slope_prev = '.')
# merge with control variables
df_ger_slope_prev <- df_ger %>% select(-time, -date, -rate_day, -socdist_single_tile) %>%
distinct() %>%
inner_join(df_ger_slope_prev, by = 'kreis') %>%
drop_na()
df_ger_slope_prev
NA
df_ger_slope_prev %>% ggplot(aes(slope_prev)) + geom_histogram(bins = 100)
df_ger_slope_socdist %>% ggplot(aes(slope_socdist)) + geom_histogram(bins = 100)
NA
NA
ctrls <- cforest_unbiased(ntree=500, mtry=5)
crf_o_fit_prev <- cforest(slope_prev ~ pers_o + women + academics +
cdu + afd + hospital_beds + gdp + manufact +
airport + age + popdens,
df_ger_slope_prev[-1],
controls = ctrls)
crf_o_varimp_prev <- varimp(crf_o_fit_prev, nperm = 1)
crf_o_varimp_cond_prev <- varimp(crf_o_fit_prev, conditional = T, nperm = 1)
ctrls <- cforest_unbiased(ntree=500, mtry=5)
crf_c_fit_prev <- cforest(slope_prev ~ pers_c + women + academics +
cdu + afd + hospital_beds + gdp + manufact +
airport + age + popdens,
df_ger_slope_prev[-1],
controls = ctrls)
crf_c_varimp_prev <- varimp(crf_c_fit_prev, nperm = 1)
crf_c_varimp_cond_prev <- varimp(crf_c_fit_prev, conditional = T, nperm = 1)
crf_c_varimp_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
ggplot(aes(x=variable, y=.)) +
geom_bar(stat = 'identity') +
theme(axis.text.x = element_text(angle = 90))
crf_c_varimp_cond_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
ggplot(aes(x=variable, y=.)) +
geom_bar(stat = 'identity') +
theme(axis.text.x = element_text(angle = 90))
ctrls <- cforest_unbiased(ntree=500, mtry=5)
crf_e_fit_prev <- cforest(slope_prev ~ pers_e + women + academics +
cdu + afd + hospital_beds + gdp + manufact +
airport + age + popdens,
df_ger_slope_prev[-1],
controls = ctrls)
crf_e_varimp_prev <- varimp(crf_e_fit_prev, nperm = 1)
crf_e_varimp_cond_prev <- varimp(crf_e_fit_prev, conditional = T, nperm = 1)
crf_e_varimp_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
ggplot(aes(x=variable, y=.)) +
geom_bar(stat = 'identity') +
theme(axis.text.x = element_text(angle = 90))
crf_e_varimp_cond_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
ggplot(aes(x=variable, y=.)) +
geom_bar(stat = 'identity') +
theme(axis.text.x = element_text(angle = 90))
ctrls <- cforest_unbiased(ntree=500, mtry=5)
crf_a_fit_prev <- cforest(slope_prev ~ pers_a + women + academics +
cdu + afd + hospital_beds + gdp + manufact +
airport + age + popdens,
df_ger_slope_prev[-1],
controls = ctrls)
crf_a_varimp_prev <- varimp(crf_a_fit_prev, nperm = 1)
crf_a_varimp_cond_prev <- varimp(crf_a_fit_prev, conditional = T, nperm = 1)
crf_a_varimp_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
ggplot(aes(x=variable, y=.)) +
geom_bar(stat = 'identity') +
theme(axis.text.x = element_text(angle = 90))
crf_a_varimp_cond_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
ggplot(aes(x=variable, y=.)) +
geom_bar(stat = 'identity') +
theme(axis.text.x = element_text(angle = 90))
ctrls <- cforest_unbiased(ntree=500, mtry=5)
crf_n_fit_prev <- cforest(slope_prev ~ pers_n + women + academics +
cdu + afd + hospital_beds + gdp + manufact +
airport + age + popdens,
df_ger_slope_prev[-1],
controls = ctrls)
crf_n_varimp_prev <- varimp(crf_n_fit_prev, nperm = 1)
crf_n_varimp_cond_prev <- varimp(crf_n_fit_prev, conditional = T, nperm = 1)
crf_n_varimp_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
ggplot(aes(x=variable, y=.)) +
geom_bar(stat = 'identity') +
theme(axis.text.x = element_text(angle = 90))
crf_n_varimp_cond_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
ggplot(aes(x=variable, y=.)) +
geom_bar(stat = 'identity') +
theme(axis.text.x = element_text(angle = 90))
social distancing ~ openness